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Abstract 

We consider the effects of noise on a model of epidemic outbreaks, where the outbreaks appear 
randomly. Using a constructive transition approach that predicts large outbreaks prior to their 
occurrence, we derive an adaptive control scheme that prevents large outbreaks from occurring. 
The theory is applicable to a wide range of stochastic processes with underlying deterministic 
structure. 
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I. INTRODUCTION 



Recently, there has been much research of steady state epidemics in random populations 
and its control j^. Non-equilibrium diseases, in contrast, are those diseases exhibiting 
outbreaks that fluctuate in time. Childhood 0, Q] and tropical diseases 0, 0] are a few ex- 
amples of outbreaks having strong annual oscillations with random amplitude. In modeling 
annual incidence of infections, random components from the environment and/or popula- 
tions play a significant role [3,0]. While excellent data from seasonally fluctuating diseases 
illustrates strong annual oscillations with random peak outbreaks in the infections jj, |^|, 
models and data analysis reveal that outbreaks stem from stochastic perturbations in either 
population or epidemic parameters, making deterministic prediction difficult. 

Predictability of seasonally driven diseases that are stochastic is necessary for the ap- 
plication of methods to sup press future outbreaks. Many vaccine schemes are available for 
equilibrium diseases P, llo| |. but in the case of non-equilibrium outbreaks, current methods 
may enhance outbreaks or fail to achieve their goals Jill 12|. (Similar problems arise in large 



fluctuation theory of stochastic dynamical systems 



14j.) Other methods pulse the popu- 



lation without sampling for prediction Q|, or rely on reducing spread via mean threshold 
reduction To address the problem of suppressing outbreaks in stochastic epidemics, we 
apply a new mathematical method to a stochastic model to predict outbreaks before 
they occur, and then adapt a vaccine strategy which prevents the outbreak from occurring. 
The theory exploits a transition probability description from small amplitude incidence to 
outbreak dynamics, and generates a region of high probability transport of the most sen- 
sitive regions to stochastic effects. Moreover, it allows us to monitor regions of stochastic 
dynamics that have a high probability of preceding a large outbreak, which in turn leads 
to a design of a vaccine control strategy to suppress outbreaks. We thus argue a general 
simple, but effective, control technique that takes advantage of complicated interactions of 
determinism and noise. The techniques introduced here may also be applied to general 
stochastic non-autonomous systems of the form 

^ = G{x,t)+r,(t) (1) 

where G(x, t) = G(x, t + 1), and the noise is added periodically with the period of drive; i.e., 

rj(t)=ri n A(t-n),n=l,2,---, (2) 



and A is the Dirac delta function, and rj n is now a discrete random variable. The form of 
Eq. [T] allows us to consider the dynamics as a discrete-time constantly perturbed stochastic 
dynamical system. 

II. A STOCHASTIC EPIDEMIC MODEL 

A standard system used to study and predict the stochastic dynamics of disease epidemics 
is based on a simplified reduced version of the well-known SEIR compartmental model 
[3, 0, Q], known as the modified SI model 1^]. In deterministic settings, the system 
has been exploited to mode, sing.e and conpled patch popu.ations Q, as we.l as testing 
vaccine strategies [13J, [19|. Assume that the population is sufficiently large so that the 
various subgroups are assumed to be continuous. The population dynamics is described by: 
Susceptible S(t). Exposed, but not yet infectious, E(t). Infective I(t). The recovered R(t) 
class in the model can be derived from model results since S + E + I + R= l 

Seasonality is input into the model via the contact rate, /3(t), so we let /3(t) = (3q(1 + 
S cos2irt), where < 8 < 1. Other parameters used to quantify the dynamics are susceptible 
input rate, fi (which includes the birth rate, as well as a possible fixed vaccine control), the 
mean latent period, a~ l , and the infectious period, 7 -1 . The full deterministic rate equations 
are given by: 

dS{t) 



dt 
dE(t) 

dt 
dl(t) 



fi(l + h(t))-P(t)SI-iiS 
P(t)SI -aE-nE (3) 
aE — 7I — /// 



dt 

R(t) = l-(S(t) + E(t)+I(t)), 

where h(t) is a small perturbation used for vaccination. That is, when h(t) is negative, the 
input of susceptibles into the system is reduced. Since it will be designed to be adaptive 
stochastic control, hit) will also depend on the state variables. 

n 

For realistic childhood disease parameters chosen here, theoretical [20] and numerical 
analysis [l7[ show that for almost all cases, the infective and exposed population follow each 
other in time to first order, leading to a reduction which describes a modified SI model 
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(MSI), given by 



dS(t) 




(4) 



dt 
dl(t) 
dt 



The parameters used for measles data [20] are given by // = 0.02, a = 1/0.0279, 7 = 1/0.01, 
Po = 1575, and 5 = 0.095, and are fixed throughout the paper. Here, the parameter h{t) is 
a time- dependent vaccine control whose value we will calculate adaptively, and depends on 
the phase space location of (S(t),I(t)). 

Following the discretized stochastic model in Eq. ^ we strobe the system with period-1 
to create a Poincare map. Without loss of generality, we define a discrete stochastic model 
for the purposes of this paper 21]. Using a discrete stochastic map approach will allow us 
to make careful and accurate interpretations in terms of the (S(t), I(t)) variables, as well 
as examine the interaction of the dynamics and control with the underlying topology of the 
system. We consider the uncontrolled stochastic system (h = 0) as a two-dimensional map, 
F, of a region D into itself 



where 77 is a two-dimensional random variable having a normal distribution given by v(x) = 



a = 0.035. Since the two dimensional deterministic system has an attractor with unequally 
sized components, the noise amplitude is scaled so that it is defined on the unit square. 
Because the standard deviation is based on the re-scaled coordinates, it is small compared to 
the attractor size and is smaller than the modulation component of the contact rate in Eq. 0] 
A typical time series of the I component is shown in Fig. ^ Notice the frequent aperiodic 
bursts, which for the chosen parameters of the deterministic part of the model, Eq. Pjl. 
would not occur were it not for the random perturbations in Eq. (JSJ); the deterministic and 
stochastic parts interact in a fundamental way to create complicated oscillations that either 
phenomenon could not create on their own. 

Notice that in the absence of any stochastic fluctuations, {f]{t) = 0), the system will 
settle down to one of two periodic solutions. The two stable solutions are plotted in the 
Fig. ^insert. The period two cycle has a small amplitude (SA) while the period 3 cycle is of 
large amplitude (LA). However, as seen from the time series in the figure, outbreaks, which 



(S,I) (t + 1) = F [(S, I) (t)]+ V (t), 



(5) 



e 



^-M/2/( 2vr ||E 



Z 2 ), with £ = diag(a 2 ), and we choose the standard deviation to be 
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occur due to stochastic fluctuations, may have enhanced amplitudes by almost an order of 
magnitude over the period 3 cycle. 

Although the system is stochastic, its dynamics may be quantified in terms of Lyapunov 



exponents by spatial integration against the invariant density [22| . For the parameters used 
to generate the time series in Fig. ^ we compute the Lyapunov exponents, and find them 
to be Ai = 0.1638 and A2 = —0.4853. These values, together with the evidence of nearly 



intersecting stable and unstable manifolds 



231 ] . indicate a comj 



eted horseshoe dynamics 
However, the corn- 




under the influence of the noise, described as stochastic chaos |8j, [1 
pleted horseshoe dynamics, indicative of chaos in deterministic systems, is a geometric way 
of thinking about the interaction of noise and the underlying manifold structure of the de- 
terministic part. The chaotic- looking dynamics is the result of mixing two stable attractors, 
while sampling unstable dynamics between them. The positive Lyapunov exponent is there- 
fore a way of measuring contributions to the stochastic attractor of dynamics tracking near 
unstable manifolds. The fraction of time spent near the unstable manifolds, as well as the 
transition probabilities of the dynamics switching from small to large amplitude behavior 
may be explained by taking a dynamic probabilistic approach, which we sketch briefly. A 
full mathematical description is given in [l5|. 



III. DISCRETE STOCHASTIC DYNAMICS AND TRANSITION PROBABILI- 
TIES 

If the noise is continuous, we can compute the evolution of the probability density using 
a Fokker-Planck approach [27| . However, since the approach is one of discrete noise as in 
Eq.^ we evolve the densities discretely as well. That is, since the solution to the periodically 
driven is computed every period to form the discrete map, we do the same with the density. 

We assume the noise comes from a distribution, u(x). The evolution of an initial probabil- 
ity density function (PDF), p : D C 3? 2 — > 3?, is defined by the stochastic Frobenius- Perron 
operator |lj P F : L 1 ^ 2 ) L 1 ^ 2 ), given by 

P F [p(x)]= [ v{x - F(y))p(y)dy. (6) 

J D 

The density is invariant if it is a fixed point of the operator. This approach allows an 
approximation of the probabilistic transitions of one part of phase space to another as 
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FIG. 1: This is an uncontrolled time series of the fraction of infectives I for the MSI model under 
random forcing in Eqs.ElandOD The parameters are given in the text. Inset: Note the small (SA) 
period 2 and large (LA) period 3 amplitude oscillations of the underlying bi-stable deterministic 
system. No chaos is present when r](t) = 0, and the system only exhibits periodic SA or LA 
oscillations. 

well as the invariant density j^J. 

To compute the transition probabilities from one region of phase space to another, we 
discretize the region D of phase space. Specifically, we assume there exists a cover of the 
region D by disjoint sets Bi, 

N 

D= \JB t . (7) 



Denning the set of characteristic basis functions, 



ipi(x) = XbA x ) = 




(8) 



allows one to generate finite dimensional projections of transport by computing the N x N 
matrix entries of a transition probability matrix Q] given by the equation 



M, 



Pf {ipi (x))ipj (x)dx. 



(9) 



D 



Therefore, Eq. El yields the probability of transporting mass from box Bi to Bj. 

In considering the problem of predicting stochastic outbreaks in the MSI model, we wish 
to compute the transition from a small amplitude (SA) oscillation to a large amplitude 
(LA) outbreak in a time series, such as the one generated in Fig. ^ The inset shows 



the deterministic periodic orbits of SA and LA, although noise may generate much larger 
outbreaks than the deterministic LA orbit. Stochastic perturbations of SA in the inset are 
approximately the same amplitude, and therefore are used as a threshold to define large 
outbreaks. The mass flux entries generated by Eq. El can be combined with the invariant 
density to generate the conditional probability of transition from set Bi to Bj, given B { . 
A representation of the transition probability is depicted in Fig. El Notice that the most 
active transport regions lie close to a stable manifold of an LA orbit (period 3 saddle) in the 
underlying deterministic system. This stable manifold is the deterministic basin boundary 
which separates the SA (period 2) and LA (period 3) regular orbits of Eq. HI and the coloring 
denotes the degree and location where this pseudo barrier is overcome due to noise. Notice 
that near each of the basin boundary saddles of period 3, transition to an outbreak is likely. 
However, the highest transition region is not near any saddle. Rather, the probability of an 
outbreak in this region is solely due to the interaction of the noise and the global topology 
of the underlying deterministic dynamics. 



IV. ACTIVE CONTROL OF STOCHASTIC OUTBREAKS 

For deterministic systems, normal methods of vaccine control will reduce the input rate of 
susceptibles. The value of h is usually computed so that at equilibrium (no seasonal forcing, 
or 5 = 0), the net rate of production of infectives in one infectious period is less than 
unity. Under these conditions, the disease will die out. However, control of small amplitude 
oscillations in the periodically driven case can be done, but the disease will persist [26]. 

In the stochastic case in presence of periodic drives, constant controls may make the 
problem worse. In Fig. El we see a direct comparison of constant vaccine control and no 
control. Notice that although the mean level of outbreaks appear to be reduced, the large 
fluctuations are greater than without control. Therefore, constant vaccine control, although 
sometimes the only guide, may increase the size of large outbreaks. Therefore, it is natural 
to try to sample and control discretely when considering stochastic outbreaks. 

Vaccine activation using a variable h depends on finding the regions where an outbreak 
is most likely upon the next iteration. These are points of the trajectory generated by Eq. El 
in the SA basin that precede iterates in the LA basin. Although we compute conditional 
outbreaks from the spatial averages using the transition matrix, this is verified temporally. 
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FIG. 2: (a) An uncontrolled time series of infective fraction as a function of time, (b) Constant 
vaccine control to reduce the rate of input of susceptibles. 



Using an uncontrolled stochastic time series of 50,000 iterates, and checking in which basin 
(SA or LA) each iterate is located, we show in Fig. 0] most likely pre-outbreak regions. In 
comparison to Fig. the spatial average predicts similar transport regions of high condi- 
tional probability of SA - LA transition. 

We now define a bull's eye (BE) region to be an open connected neighborhood having high 
probability of transition from SA to LA outbreaks. The BE region, for a chosen threshold, 
is clearly shown in red in Fig. El Distinguishing the center point x c , the BE region includes 
a neighborhood of radius e that has a probability greater than a given threshold. Notice 
that this is not the only region in which transition occurs. Monitoring the BE region alone, 
therefore, is not sufficient for prediction of transitions [25]. However, it can be used to 
determine other regions that are not obvious for transition to an outbreak. 

We can use the BE region as a first guess to monitor the dynamics. Let xq G E be the 
current point of the observed dynamics, and xl a desired target point in the transport space 
close to the image of E, but in a region of lower transition probability. The relationship 
between the current point in the trajectory x and the center of the bull's eye x c is x = x c +y 
for some y. To move the image of xq closer to the target point xl, we activate the control 
parameter h in Eq. |3J By Taylor expansion about F(x c , fi) when h = and ignoring higher 



S 



FIG. 3: The GTM result of the conditional probability of transition from small amplitudes to large 
outbreaks using the same parameters as in Fig. ^ The highest probability regions of transport 
(red) point to a bull's eye monitoring region for control. Overlaid are the stable and unstable 
manifolds corresponding to the underlying deterministic model. 
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FIG. 4: (Color) Temporal average of those iterates leading to outbreaks in the next iterate using 
the same parameters as in Fig. ^ Notice the agreement with spatial average in Fig. |31 The ellipse 
bounds the detection region. The figure eight curve is the image of the ellipse with controlled 
targeting. See text for details. 

order terms, we solve, 

_ - F(x c , ji) - d x F(x c , /i)y) T ^F(x c , fi) 

\\d,F(x c ^W [ } 

assuming d IJi F(x c , fi) ^ 0. This control strategy is designed to target a desired region of 




lower probability, given the iterates land in a region of high transition probability. 

Now we apply control to suppress large amplitude outbreaks. Focusing on points in a 
neighborhood of the BE region has the disadvantage that the values of I are already fairly 
large. Therefore, we use the detection region of the neighborhood around the (deterministic) 
pre-image of the bull's eye, F~ 1 (E, /x), shown as an ellipse in Fig.|3J Using Eq. the image 
of the ellipse, Ih, is found to be the figure eight shown in Fig. HJ We targeted a region in 
Ih which is close to the BE region but has a very low transition probability. Our techniques 
successfully steer trajectories away from the bull's eye region towards SA behavior by using 
only vaccine perturbations that control the flow of susceptibles about some mean value. 

One advantage of choosing the detection region to be the pre-image of BE is for relatively 
low values for the number of infected individuals (J), a prediction can be made about the 
future increase and steps can be taken to avert these dynamics. The perturbations represent 
a vaccination program, taking the form of \x new = + h(t)). If h is negative, then more 
vaccinations are required to reduce the rate of susceptible individuals being introduced 
into the population. An example of the success of this algorithm is shown in Fig. ED On 
average, perturbations are applied 25-30% of the time. Notice the maximum amplitude 
in comparison to the uncontrolled dynamics of Fig. For this example, the Lyapunov 
exponents are Ai = 0.0794 and A2 = —0.3764, where the maximum exponent has been 
significantly decreased. 

V. DISCUSSION 

Stochastic bursting is present in many systems that are based on population dynamic 
modeling. In general, when such systems are subject to periodic forcing, there exist param- 
eter regions in which multiple attractors co-exist. Typically, one of these attractors arises 
from periodically forced equilibrium, and therefore, is typically of small amplitude. On the 
other hand, the other attractors bifurcate from saddle node orbits, which tend be of larger 
amplitude. Such bi-stable systems can have simple manifold structure, but when considered 
in the presence of stochastic fluctuations, may exhibit complex mixing between the bi-stable 
attractors, coupled with complicated looking transients between the basins. 

By using the PDF flux, we are able to distinguish regions in the small amplitude basin 
that are quite sensitive to stochastic effects. We use this information in a control algorithm 
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FIG. 5: Stochastic control to suppress large outbreaks in the MSI model, (a) Infectives with 
suppressed outbreaks due to control in the influx of infectives, (b) Perturbations h to the susceptible 
input rate \x in Eq. 0J 

to prevent bursting dynamics (that is, to control stochastic chaos). It monitors this sensitive 
region and adjusts one physically relevant parameter to keep trajectories in the SA basin. 
This idea of monitoring a loss region has been used in other chaos control schemes that 



are deterministic (i.e. 



28, 



3). 



To our knowledge, we are not aware of any stochastic chaos 



control methods that account specifically for the emergent effects of stochastic perturbations. 

One concern with a probabilistic detection scheme is that it is dependent on the choice 
of monitoring region used for transition to an outbreak. Two issues with taking an actual 
time series and using the monitoring scheme above is that it may miss an outbreak that is 
there (missed detection), or it may predict an outbreak that does not occur. These statistics 
depend heavily on the size of the monitoring one uses. To see this in Fig. |Ul we change the 
radius around of the center of the bull's eye and the radius around it's preimage. Each dot 
plotted in Fig. |Hl is for a different radius. The smallest radii are represented by the data 
points on the right. As we increase the radii, the data points move along the curve to the 
left. The false alarms are those outbreaks predicted by the bull's eye, but do not occur. 
The missed detection are the bursts that occur but are not predicted by the maximum flux 
hypothesis. It is the percentage not detected. 

The choice we made for the detection region has solely been guided by time series ob- 
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FIG. 6: A plot showing optimized predictability. Plotted are the false alarm rates versus missed 
detections. See text for details. 

servations and PDF flux predictions. It has not been optimized for the minimum number 
or size of perturbations. Because of the stochastic perturbations added to the system, the 
control measures will not "trail off" as in targeting unstable periodic orbits. 

One very interesting aspect of the control perturbation is that the value of h is sometimes 
positive. This counter-intuitive result may be explained since an increase in S is used to 
trigger an earlier, but smaller outbreak. To understand this, we consider the MSI model 
but transformed and scaled, so that the steady state equilibrium in the absence of forcing 
is now at the origin, and we examine the conservative system in the absence of damping as 
well 0: 

x\t) = -vy (11) 
y'(t) = vx(l + y). 

In Eq. HU x is a scaled susceptible, y is a scaled infective, the equilibrium is at the origin, 
and the frequency v is a function of the epidemiological parameters. Notice that since the 
population is assumed to be constant, in the absence of any infectives (y — — 1 in Eq. ITTj) . 
the fraction of susceptibles slowly increases. In addition, all oscillatory solutions must lie on 
level curves to the Lyapunov function: V(x, y) = x 2 + 2y — 2ln(y + 1). 

Now suppose we have a small amount of infectives imposed by a strong level of vaccine. 
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Then the infectives will stay small for a long period of time, until enough susceptibles grow 
to cause an outbreak of very large amplitude by coming in contact with a few infectives [31]. 
That is, an outbreak will not occur unless the susceptibles reach a critical level in a long 
time scale while in the presence of a small fraction of infectives. To be specific, suppose 
y = — 1 + ce, where c > is constant. Then y' = xce. If x < 0, then the infectives decrease 
further, implying a much larger outbreak at a later time. Therefore, if one increases the 
infectives, the system fires sooner, with a smaller outbreak, since the infectives are pushed 
further away from the invariant line y — — 1. When the control, h, is adjusted so that it is 
positive, the effect is to cause an increase in the rate of infectives, thus reducing the size of 
the outbreak. 

Finally, although the vaccine control fluctuations do not decrease the mean incidence 
evels of infection, the control may be combined with tracking methods for epidemic control 



26| to reduce the mean reproductive rate of infection below threshold to kill off the disease 



without causing unwanted outbreaks during vaccination. 
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